rm(list = ls(all = TRUE))
gc()
##           used (Mb) gc trigger (Mb) max used (Mb)
## Ncells  555850 29.7    1244914 66.5   686457 36.7
## Vcells 1025182  7.9    8388608 64.0  1877160 14.4

1 packages

# %>% 
library(magrittr)
library(ggplot2)


`%!in%` = Negate(`%in%`)

2 palette

n = 4
col = RColorBrewer::brewer.pal(n, 'Paired')

pal = col[seq(1, n, 2)]
my_dir = file.path("..", "reports", "Table7")
if (!dir.exists(my_dir)) {
  dir.create(my_dir)
}
fr = file.path('..', 'output',  'Table7_stat.txt')
file.create(fr)
## [1] TRUE

3 f-ctions

3.1 summary stat

summary_stats <- function(df, measurevar, groupvars, conf.level = 0.95) {
  df %>%
    dplyr::group_by(across(all_of(groupvars))) %>%
    dplyr::summarise(
      N = sum(!is.na(.data[[measurevar]])),
      mean = mean(.data[[measurevar]], na.rm = TRUE),
      sd = sd(.data[[measurevar]], na.rm = TRUE),
      .groups = "drop"
    ) %>%
    dplyr::mutate(
      se = sd / sqrt(N),
      ci = se * qt(conf.level/2 + 0.5, N - 1)
    )
}

3.2 scale to the average of the reference group

process_data <- function(df, groupvars, measurevar, scale_treatment = "noninoculated") {

  df.long = df %>%
    tidyr::pivot_longer(cols = 3:ncol(df), names_to = "transcript",
                        values_to = measurevar, values_drop_na = TRUE)
  data.SE = summary_stats(df.long, measurevar, groupvars)
  scale_reference = data.SE %>%
    dplyr::filter(Treatment == scale_treatment) %>%
    dplyr::select(Tissue, transcript, mean) %>%
    dplyr::rename(scale_mean = mean)
  df.scaled <- dplyr::left_join(df.long, scale_reference, by = c("Tissue", "transcript")) %>%
    dplyr::mutate(scaled = .data[[measurevar]] / scale_mean) %>%
    dplyr::select(-scale_mean)
  df.scaled = df.scaled %>% 
    dplyr::arrange(dplyr::desc(transcript), dplyr::desc(Treatment))
  shoots = df.scaled %>% dplyr::filter(Tissue == "shoots")
  roots = df.scaled %>% dplyr::filter(Tissue == "roots")
  
  return(list(
    shoots = shoots,
    roots = roots
  ))
}

3.3 permutation t-test

  • for non-i.i.d, hypotheses
  • a non-parametric robust alternative to traditional t-tests
  • small sample sizes
  • robust to unequal variances
perm_test_by_transcript <- function(mydata.long, measurevar = "measurement", groupvar = "Treatment") {
  
  temp = data.frame()
  transcript_levels = levels(mydata.long$transcript)
  treatment_levels = levels(mydata.long[[groupvar]])
  treatment_pairs = combn(treatment_levels, 2, simplify = FALSE)
  
  for(i in transcript_levels) {
    data = mydata.long[mydata.long$transcript == i, ]
    k = 8 # nrow(data); an exact permutation test, assuming the groups are balanced (i.e., 4 observations in each group)
    pvalues = purrr::map_dbl(treatment_pairs, function(x) {
      subset_data = base::subset(data, data[[groupvar]] %in% x)
      res = MKinfer::perm.t.test(
        formula = stats::as.formula(base::paste(measurevar, "~", groupvar)),
        data = subset_data,
        alternative = "two.sided",
        mu = 0,
        paired = FALSE,
        var.equal = FALSE,
        conf.level = 0.95,
        perm.conf.int = 0.95,
        symmetric = TRUE,
        p.adjust.method = "holm",
        detailed = TRUE,
        R = choose(k, k/2), # sum(choose(k, 1:(k-1))) for unbalanced groups
        set.seed = 123456
      )
      res$perm.p.value
    })
    
    tmp = base::as.data.frame(base::t(pvalues))
    colnames(tmp) = purrr::map_chr(treatment_pairs, ~base::paste(.x, collapse = ' vs '))
    rownames(tmp) = i
    
    temp <- base::rbind(temp, tmp)
  }
  
  return(temp)
}

3.4 plot

plot_gene_expression <- function(data.SE
                                 , data.long
                                 , stat.test.sig
                                 , transcripts_excl = c('13-LOX', 'PTI5')
                                 , color_values
                                 , facet_cols = 6
                                 , y_scales = NULL
                                 , dodge_width = 0.8,
                                 plot_title = ""
                                 ) {
  
  dodge = position_dodge(width = dodge_width)
  
  data.SE.filtered = dplyr::filter(data.SE, !transcript %in% transcripts_excl)
  data.long.filtered = dplyr::filter(data.long, !transcript %in% transcripts_excl)
  
  stat.test.filtered = NULL
  if (nrow(stat.test.sig) > 0) stat.test.filtered = dplyr::filter(stat.test.sig, !transcript %in% transcripts_excl)
  
  p = ggplot(data.SE.filtered, aes(x = Treatment, y = mean)) +
    geom_point(size=3.5, shape = 22, position = dodge, aes(fill = Treatment), colour = "black") +
    geom_point(data = data.long.filtered, aes(x = Treatment, y = scaled, fill = Treatment), colour = "black",
               size = 2.0, shape = 21,
               position = dodge) +
    geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.3, lwd = 0.5,
                  position = dodge) +
    facet_wrap(~ transcript, ncol = facet_cols, scales = "free", drop = TRUE)
  
  if (!is.null(y_scales)) {
    p = p + ggh4x::facetted_pos_scales(y = y_scales)
  }
  
  p = p +
    # geom_hline(yintercept = 2.5, alpha = 0.0) +
    geom_hline(yintercept = 0, alpha = 0.0) +
    geom_hline(yintercept = 1, alpha = 0.5, linetype = "dotted", col = "gray45") +
    ggtitle(plot_title) +
    theme_bw() +
    scale_colour_manual(name = "", values = rev(color_values)) +
    scale_fill_manual(name = "", values = rev(color_values)) +
    labs(x = "", y = "Relative gene expression (+/- SE)"
         # , subtitle = "Permutation t-test measurements"
         ) +
    theme(plot.subtitle = element_text(size = 10),
          axis.text = element_text(size = 12.5),
          axis.text.x = element_text(size = 12.5, angle = 90),
          axis.title = element_text(size = 12.5, face = "bold"),
          strip.text = element_text(size = 12.5),
          title = element_text(size = 15, face = "bold"),
          # axis.ticks.x = element_blank(),
          legend.key.height = unit(1.5, "cm"),
          legend.key.width = unit(1.75, "cm"),
          legend.text = element_text(size = 12.5),
          legend.title = element_text(size = 12.5),
          legend.background = element_rect(fill = "transparent", size = 0.5, linetype = "dotted"),
          legend.position = "top",
          legend.justification = "right",
          plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.title.y.right = element_blank(),
          axis.text.y.right = element_blank(),
          # axis.ticks.y = element_blank(),
          axis.text.y = element_text(size = 12.5, margin = margin(r = 0)),
          panel.spacing = unit(1, "lines"),
          strip.background = element_rect(size = 0.5, fill = "transparent", color = NA) ,
          panel.border = element_blank(),
          axis.line = element_line(color = "black")
          )
  
    p = p + theme(legend.position = "none")

  


  
  if (!is.null(stat.test.filtered) && nrow(stat.test.filtered) > 0) {
    
    p = p + ggpubr::stat_pvalue_manual(
      stat.test.filtered,
      label = "p.adj.signif",
      xmin = "xmin",
      xmax = "xmax",
      y.position = "y.position",
      hide.ns = FALSE,
      tip.length = 0.01,
      step.increase = 0.0,
      inherit.aes = FALSE 
    )
    
  }
    
  return(p)
  
}

3.5 distribution, variance and effects

  • on qPCR data
  • Denote: LOQ values skew distributions
  • Denote: distribution tests: low power for small sample size, see QQ of residuals
  • if most or all points fall inside the shaded confidence band, the sample’s distribution does not show strong evidence of departure from normality at the plotted sample size and confidence level
  • a few isolated points outside the band at the extremes are common with small samples and do not necessarily indicate a severe problem
  • Denote: heteroscedasticity: some tests are meant to be used with normally distributed data, but can tolerate relatively low deviation from normality; Levene’s test with mean, more robust test Brown-Forsythe with median, Fligner-Killeen when data are non-normally distributed or when problems related to outliers in the dataset cannot be resolved
  • Denote: effect size: Cohen’s d - parametric with assumptions, Wilcoxon Effect Size - the non-parametric alternative; different sensitivity to outliers
dens_and_effect <- function(mydata.long, p_colors) {
  
  # Plot density with x="measurement"
  p1 = ggpubr::ggdensity(mydata.long, x = "measurement",
              add = "mean", rug = TRUE,
              color = "Treatment", fill = "Treatment",
              facet.by = 'transcript') +
    scale_fill_manual(name = "", values = rev(p_colors)) +
    scale_color_manual(name = "", values = rev(p_colors)) +
    ggtitle("measurement") +
    facet_wrap(~transcript, ncol = 4, scales = "free")
  print(p1)
  
  # optional: Plot density with x="scaled"- sam e shape wih a shift
  # p2 = ggpubr::ggdensity(mydata.long, x = "scaled",
  #             add = "mean", rug = TRUE,
  #             color = "Treatment", fill = "Treatment",
  #             facet.by = 'transcript') +
  #   scale_fill_manual(name = "", values = rev(p_colors)) +
  #   scale_color_manual(name = "", values = rev(p_colors)) +
  #   ggtitle("scaled") +
  #   facet_wrap(~transcript, ncol = 4, scales = "free")
  # print(p2)
  
  
  cat(crayon::red('####  ####  \nDistribution tests\n####  ####  \n'))
  
    distr_tests = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      x = na.omit(.x$measurement)
      if (length(x) < 3) {
        return(tibble::tibble(
          Shapiro_Wilk = NA_real_,
          Anderson_Darling = NA_real_,
          Lilliefors_KS = NA_real_,
          Jarque_Bera = NA_real_,
          DAgostino_Skewness = NA_real_,
          n = length(x),
          transcript = .y$transcript
        ))
      }
      tibble::tibble(
        Shapiro_Wilk = tryCatch(shapiro.test(x)$p.value, error = function(e) NA_real_),
        Anderson_Darling = tryCatch(nortest::ad.test(x)$p.value, error = function(e) NA_real_),
        Lilliefors_KS = tryCatch(nortest::lillie.test(x)$p.value, error = function(e) NA_real_),
        Jarque_Bera = tryCatch(tseries::jarque.bera.test(x)$p.value, error = function(e) NA_real_),
        DAgostino_Skewness = tryCatch(moments::agostino.test(x)$p.value, error = function(e) NA_real_),
        n = length(x),
        transcript = .y$transcript
      )
    }) %>%
    dplyr::bind_rows() %>%
    dplyr::ungroup() %>%
    dplyr::mutate(
      Shapiro_Wilk_BH = stats::p.adjust(Shapiro_Wilk, method = "BH"),
      Anderson_Darling_BH = stats::p.adjust(Anderson_Darling, method = "BH"),
      Lilliefors_KS_BH = stats::p.adjust(Lilliefors_KS, method = "BH"),
      Jarque_Bera_BH = stats::p.adjust(Jarque_Bera, method = "BH"),
      DAgostino_Skewness_BH = stats::p.adjust(DAgostino_Skewness, method = "BH")
    )

  

  combined_results = distr_tests[grep("^n$|transcript|_BH", colnames(distr_tests))]
  print(combined_results)

  
  

  cat(crayon::red("####  ####  \nQuantile-Quantile plots\n####  ####  \n"))
  qq_plot_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      ggpubr::ggqqplot(.x$measurement) +
        ggplot2::ggtitle(paste("Raw:", .y$transcript))
    })
  
  resid_qq_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      model = lm(measurement ~ Treatment, data = .x)
      resids = resid(model)
      ggpubr::ggqqplot(resids) +
        ggplot2::ggtitle(paste("Residuals:", .y$transcript))
    })
  

  resid_density_list = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_map(~ {
      model = lm(measurement ~ Treatment, data = .x)
      resids = scale(resid(model))  # standardize residuals
      df = data.frame(resids = as.numeric(resids))
      ggplot(df, aes(x = resids)) +
        geom_density(fill = "steelblue", alpha = 0.6) +
        stat_function(fun = dnorm, color = "red", linetype = "dashed") +
        ggtitle(paste("Residuals:", .y$transcript)) +
        theme_minimal()
    })
  
  # Arrange and print all plots
  qq_arranged = ggpubr::ggarrange(plotlist = qq_plot_list, ncol = 4, nrow = ceiling(length(qq_plot_list) / 4))
  resid_qq_arranged = ggpubr::ggarrange(plotlist = resid_qq_list, ncol = 4, nrow = ceiling(length(resid_qq_list) / 4))
  resid_hist_arranged = ggpubr::ggarrange(plotlist = resid_density_list, ncol = 4, nrow = ceiling(length(resid_density_list) / 4))
  
  print(qq_arranged)
  print(resid_qq_arranged)
  print(resid_hist_arranged)
  
  # Return invisibly
  invisible(list(qq = qq_arranged, resid_qq = resid_qq_arranged, resid_hist = resid_hist_arranged))



  
  cat(crayon::red("####  ####  \nTest for homogeneity of variance across groups\n####  ####  \n"))
  
  cat(crayon::red("####  ####  \nLevene\n####  ####  \n"))
  # Levene's test on measurement
  lev_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::levene_test(measurement ~ Treatment)
  print(lev_meas)
  
  cat(crayon::red("####  ####  \nBrown-Forsythe\n####  ####  \n"))
  # center = "median" switches Levene’s test to the Brown-Forsythe variant
  bf_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::levene_test(measurement ~ Treatment, center = "median")
  print(bf_meas)
  
  cat(crayon::red("####  ####  \nFligner\n####  ####  \n"))
  fk_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    dplyr::group_modify(~ broom::tidy(stats::fligner.test(measurement ~ Treatment, data = .x)))
  print(fk_meas)


  

  cat(crayon::red("####  ####  \nWilcoxon effect size\n####  ####  \n"))
  
  # Wilcoxon effect size on measurement
  eff_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::wilcox_effsize(measurement ~ Treatment)
  print(eff_meas)
  
  # Wilcoxon effect size on scaled
  # eff_scaled = mydata.long %>%
  #   dplyr::group_by(transcript) %>%
  #   rstatix::wilcox_effsize(scaled ~ Treatment)
  # print(eff_scaled)
  
  cat(crayon::red("####  ####  \nCohen's d Measure of Effect Size\n####  ####  \n"))
  
  # Cohen's d on measurement
  coh_meas = mydata.long %>%
    dplyr::group_by(transcript) %>%
    rstatix::cohens_d(measurement ~ Treatment, paired = FALSE)
  print(coh_meas)
  

  invisible(list(
    p_measurement = p1,
    distr_measurement = combined_results,
    levene_measurement = lev_meas,
    brownforsythe_measurement = bf_meas,
    fligner_measurement = fk_meas,
    wilcoxon_measurement = eff_meas,
    cohensd_measurement = coh_meas
  ))
}

3.6 ggplot limits

make_scale_params <- function(maxval) {
  upper = ifelse(maxval > 2, maxval + 2, maxval + 0.1) # 0.5
  step = dplyr::case_when(
    maxval >= 50 ~ 25,
    maxval >= 30 ~ 10,
    maxval >= 15 ~ 5,
    maxval >= 10 ~ 2,
    maxval >  2  ~ 1,
    TRUE         ~ 0.5
  )
  upper_round = ceiling(upper / step) * step
  breaks = seq(0, upper_round, by = step)
  list(limits = c(0, upper_round), breaks = breaks)
}

build_y_scales_for <- function(names_vec, max_per_transcript) {
  max_per_transcript %>%
    dplyr::filter(transcript %in% names_vec) %>%
    purrr::pmap(function(transcript, max_scaled) {
      params = make_scale_params(max_scaled)
      lhs = rlang::expr(transcript == !!transcript)
      rhs = rlang::expr(ggplot2::scale_y_continuous(limits = !!params$limits, breaks = !!params$breaks))
      rlang::new_formula(lhs, rhs)
    })
}

3.7 test and plot

test_and_plot <- function(data_long_raw
                          , myorder
                          , pal
                          , what
                          , plot_gene_expression_func
                          , groupvars = c("Treatment", "transcript")
                          , y_scales1
                          , y_scales2
                          ) {

  
  mydata.long = dplyr::as_tibble(data.table::data.table(data_long_raw))
  mydata.long$transcript = factor(mydata.long$transcript, levels = myorder)
  mydata.long = mydata.long %>% dplyr::arrange(factor(transcript, levels = myorder))
  
  mydata.long.SE = summary_stats(mydata.long, measurevar = "scaled", groupvars = groupvars)
  
  results = dens_and_effect(mydata.long, p_colors = pal)
  
  cat(what, file = fr, append = TRUE, sep = "\n")
  
  cat("\nDistribution tests", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$distr_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$distr_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  
  cat("\nLevene's test for homogeneity of variance across groups", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$levene_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$levene_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("\nBrown-Forsythe (Levene with median)", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$brownforsythe_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$brownforsythe_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("\nFligner-Killeen test of homogeneity of variances", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$fligner_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$fligner_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  
  cat("\nWilcoxon effect size", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$wilcoxon_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$wilcoxon_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("\nCohen's d Measure of Effect Size", file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(results$cohensd_measurement), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(results$cohensd_measurement, 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  
  cat("", file = fr, append = TRUE, sep = "\n")


  
  temp = perm_test_by_transcript(mydata.long, measurevar = "measurement", groupvar = "Treatment")
  temp$transcript = rownames(temp)
  perm = tidyr::gather(temp, contrast, perm.p, colnames(temp)[1]:colnames(temp)[ncol(temp)-1], factor_key=TRUE)
  perm$group1 = gsub(' vs.*', '', perm$contrast)
  perm$group2 = gsub('.* vs ', '', perm$contrast)
  perm$perm.p.adj = p.adjust(perm$perm.p, method = 'BH')
  perm$perm.p.adj.signif = 'ns'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.0001] = '**'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.001] = '**'
  perm$perm.p.adj.signif[perm$perm.p.adj < 0.05] = '*'
  
  
  stat.test = perm %>%
    dplyr::select(transcript, group1, group2, perm.p, perm.p.adj, perm.p.adj.signif) %>%
    dplyr::rename(
      p = perm.p,
      p.adj = perm.p.adj,
      p.adj.signif = perm.p.adj.signif
    )
  
  # Compute y-position using both mean ± se and scaled values
  y_pos_df = mydata.long.SE %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(
      max_mean_se = max(mean + se, na.rm = TRUE),
      min_mean_se = min(mean - se, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    dplyr::left_join(
      mydata.long %>%
        dplyr::group_by(transcript) %>%
        dplyr::summarise(
          max_scaled = max(scaled, na.rm = TRUE),
          min_scaled = min(scaled, na.rm = TRUE),
          .groups = "drop"
        ),
      by = "transcript"
    ) %>%
    dplyr::mutate(
      plot_max = pmax(max_mean_se, max_scaled, na.rm = TRUE),
      plot_min = pmin(min_mean_se, min_scaled, na.rm = TRUE),
      plot_range = pmax(plot_max - plot_min, 1e-6)
    ) %>%
    dplyr::select(transcript, plot_max, plot_range)
  
  # Map group names to numeric x positions
  pos_map = mydata.long %>%
    dplyr::distinct(transcript, Treatment) %>%
    dplyr::mutate(xpos = as.numeric(factor(Treatment, levels = levels(mydata.long$Treatment))))
  
  # Attach y-position and x-position info to stat.test
  stat.test = stat.test %>%
    dplyr::left_join(y_pos_df, by = "transcript") %>%
    dplyr::left_join(pos_map %>% dplyr::rename(group1 = Treatment, xmin = xpos), by = c("transcript", "group1")) %>%
    dplyr::left_join(pos_map %>% dplyr::rename(group2 = Treatment, xmax = xpos), by = c("transcript", "group2")) %>%
    dplyr::mutate(
      xmin = dplyr::if_else(is.na(xmin), as.numeric(factor(group1, levels = levels(mydata.long$Treatment))), xmin),
      xmax = dplyr::if_else(is.na(xmax), as.numeric(factor(group2, levels = levels(mydata.long$Treatment))), xmax)
    )
  
  # Compute distinct y.position per comparison for significant results
  stat.test.sig = stat.test %>%
    dplyr::filter(!is.na(p.adj) & p.adj <= 0.05) %>%
    dplyr::group_by(transcript) %>%
    dplyr::arrange(xmin, xmax, .by_group = TRUE) %>%
    dplyr::mutate(
      base_y = plot_max + 0.03 * plot_range,
      inc = 0.05 * plot_range,
      y.position = base_y + (dplyr::row_number() - 1) * inc
    ) %>%
    dplyr::ungroup() %>%
    dplyr::select(transcript, group1, group2, p.adj, p.adj.signif, y.position, xmin, xmax)

  
  group1 =  sapply(y_scales1, function(x) {
    s = deparse(x)
    s_full = paste(s, collapse = " ")
    sub('.*transcript == *"([^"]+)".*', '\\1', s_full)
  })
  group2 =  sapply(y_scales2, function(x) {
    s = deparse(x)
    s_full = paste(s, collapse = " ")
    sub('.*transcript == *"([^"]+)".*', '\\1', s_full)
  })


  p1 = plot_gene_expression_func(
    data.SE = mydata.long.SE,
    data.long = mydata.long,
    stat.test.sig = stat.test.sig,
    transcripts_excl = group2,
    facet_cols = 2,
    color_values = pal,
    plot_title = what,
    y_scales = y_scales1,
    dodge_width = 0.8
  )
  print(p1)
  
  p2 = plot_gene_expression_func(
    data.SE = mydata.long.SE,
    data.long = mydata.long,
    stat.test.sig = stat.test.sig,
    transcripts_excl = group1,
    facet_cols = 6,
    color_values = pal,
    plot_title = what,
    y_scales = y_scales2,
    dodge_width = 0.8
  )
  print(p2)
  
  return(list(plot1 = p1, plot2 = p2, stat.test = stat.test))
}

4 Test

4.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "AMFcolonisation.xlsx" "qPCRmycorrhiza.xlsx"  "README.md"           
## [4] "Table S1.xlsx"        "Table S2.xlsx"        "Table S7.xlsx"
fn = 'Table S7.xlsx'

PS216 = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = 'Test',
                            startRow = 1,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

data.table::setDT(PS216)
PS216[, Sample.ID := NULL]
PS216[, Genotype := NULL]
PS216$Tissue = as.factor(trimws(PS216$Tissue))
PS216$Treatment = factor(trimws(PS216$Treatment), levels = c("noninoculated", "inoculated"))

PS218 = PS216[grep('PS-216', PS216$Strain, invert = TRUE), ]
PS216 = PS216[grep('PS-218', PS216$Strain, invert = TRUE), ]
PS216[, Strain := NULL]
PS218[, Strain := NULL]

4.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("RBOHD", "PR1B", "CPI8", "CAB", "BGLU2", "HSP70", "PTI5", "13-LOX")
group1_names = c("PTI5", "13-LOX")

4.3 results

shoots.216 = process_data(PS216, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
roots.216 = process_data(PS216, groupvars, measurevar, scale_treatment = "noninoculated")$roots
shoots.218 = process_data(PS218, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
roots.218 = process_data(PS218, groupvars, measurevar, scale_treatment = "noninoculated")$roots


cat("Test", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")

4.3.1 shoots 216

max_per_transcript = shoots.216 %>%
  dplyr::group_by(transcript) %>%
  dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
group2_names = setdiff(max_per_transcript$transcript, group1_names)
y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
y_scales2 = build_y_scales_for(group2_names, max_per_transcript)

result = test_and_plot(data_long_raw = shoots.216
                       , myorder = myorder
                       , pal = pal
                       , what = "shoots 216"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## ####  ####  
## Distribution tests
## ####  ####
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 RBOHD              0.829               0.749            0.804  
## 2     8 PR1B               0.00129             0.00125          0.00813
## 3     8 CPI8               0.829               0.749            0.804  
## 4     8 CAB                0.829               0.749            0.804  
## 5     8 BGLU2              0.829               0.749            0.804  
## 6     8 HSP70              0.829               0.749            0.804  
## 7     8 PTI5               0.495               0.427            0.264  
## 8     8 13-LOX             0.0905              0.108            0.131  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6   0.181   0.685
## 2 PR1B           1     6   0.478   0.515
## 3 CPI8           1     6   0.725   0.427
## 4 CAB            1     6   0.00927 0.926
## 5 BGLU2          1     6   0.739   0.423
## 6 HSP70          1     6   0.0138  0.910
## 7 PTI5           1     6   1.69    0.242
## 8 13-LOX         1     6   2.61    0.157
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6   0.181   0.685
## 2 PR1B           1     6   0.478   0.515
## 3 CPI8           1     6   0.725   0.427
## 4 CAB            1     6   0.00927 0.926
## 5 BGLU2          1     6   0.739   0.423
## 6 HSP70          1     6   0.0138  0.910
## 7 PTI5           1     6   1.69    0.242
## 8 13-LOX         1     6   2.61    0.157
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 RBOHD         0.0955   0.757         1 Fligner-Killeen test of homogeneity of…
## 2 PR1B          0.0945   0.759         1 Fligner-Killeen test of homogeneity of…
## 3 CPI8          0.120    0.729         1 Fligner-Killeen test of homogeneity of…
## 4 CAB           0.0955   0.757         1 Fligner-Killeen test of homogeneity of…
## 5 BGLU2         0.123    0.726         1 Fligner-Killeen test of homogeneity of…
## 6 HSP70         0.0937   0.760         1 Fligner-Killeen test of homogeneity of…
## 7 PTI5          1.45     0.228         1 Fligner-Killeen test of homogeneity of…
## 8 13-LOX        1.46     0.227         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.102 RBOHD          4     4 small    
## 2 measurement noninoculated inoculated   0.306 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated   0.510 CPI8           4     4 large    
## 4 measurement noninoculated inoculated   0.102 CAB            4     4 small    
## 5 measurement noninoculated inoculated   0.612 BGLU2          4     4 large    
## 6 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.714 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.269 RBOHD          4     4 small    
## 2 measurement noninoculated inoculated  -0.696 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated  -1.21  CPI8           4     4 large    
## 4 measurement noninoculated inoculated  -0.425 CAB            4     4 small    
## 5 measurement noninoculated inoculated  -1.44  BGLU2          4     4 large    
## 6 measurement noninoculated inoculated  -0.389 HSP70          4     4 small    
## 7 measurement noninoculated inoculated  -4.45  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -1.68  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result$stat.test))]
print(res[res$p.adj.signif != 'ns', ])
## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)
cat("", file = fr, append = TRUE, sep = "\n")
output_text = "permutational t-test"
cat(output_text, file = fr, append = TRUE, sep = "\n")
header = paste(colnames(res), collapse = "\t")
cat(header, file = fr, append = TRUE, sep = "\n")
output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")

  

ggsave(
  filename = file.path(my_dir, "Test_shoots.216_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_shoots.216_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

4.3.2 roots 216

max_per_transcript = roots.216 %>%
  dplyr::group_by(transcript) %>%
  dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
group2_names = setdiff(max_per_transcript$transcript, group1_names)
y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
y_scales2 = build_y_scales_for(group2_names, max_per_transcript)

result = test_and_plot(data_long_raw = roots.216
                       , myorder = myorder
                       , pal = pal
                       , what = "roots 216"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 RBOHD             0.500               0.532            0.619   
## 2     8 PR1B              0.000177            0.000309         0.000230
## 3     8 CPI8              0.0915              0.121            0.258   
## 4     8 CAB               0.0543              0.0534           0.0216  
## 5     8 BGLU2             0.000188            0.000339         0.000616
## 6     8 HSP70             0.765               0.809            0.956   
## 7     8 PTI5              0.0915              0.126            0.258   
## 8     8 13-LOX            0.500               0.532            0.580   
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 RBOHD          1     6    1.19   0.317 
## 2 PR1B           1     6    1.00   0.356 
## 3 CPI8           1     6    7.39   0.0347
## 4 CAB            1     6    3.90   0.0956
## 5 BGLU2          1     6    0.792  0.408 
## 6 HSP70          1     6    0.0123 0.915 
## 7 PTI5           1     6    9.03   0.0238
## 8 13-LOX         1     6    0.258  0.630 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 RBOHD          1     6    1.19   0.317 
## 2 PR1B           1     6    1.00   0.356 
## 3 CPI8           1     6    7.39   0.0347
## 4 CAB            1     6    3.90   0.0956
## 5 BGLU2          1     6    0.792  0.408 
## 6 HSP70          1     6    0.0123 0.915 
## 7 PTI5           1     6    9.03   0.0238
## 8 13-LOX         1     6    0.258  0.630 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 RBOHD       1.47      0.226          1 Fligner-Killeen test of homogeneity of…
## 2 PR1B        4.99      0.0255         1 Fligner-Killeen test of homogeneity of…
## 3 CPI8        2.85      0.0913         1 Fligner-Killeen test of homogeneity of…
## 4 CAB         5.00      0.0253         1 Fligner-Killeen test of homogeneity of…
## 5 BGLU2       0.619     0.431          1 Fligner-Killeen test of homogeneity of…
## 6 HSP70       0.624     0.430          1 Fligner-Killeen test of homogeneity of…
## 7 PTI5        4.96      0.0259         1 Fligner-Killeen test of homogeneity of…
## 8 13-LOX      0.000202  0.989          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.306 RBOHD          4     4 moderate 
## 2 measurement noninoculated inoculated   0.102 PR1B           4     4 small    
## 3 measurement noninoculated inoculated   0.408 CPI8           4     4 moderate 
## 4 measurement noninoculated inoculated   0.408 CAB            4     4 moderate 
## 5 measurement noninoculated inoculated   0.306 BGLU2          4     4 moderate 
## 6 measurement noninoculated inoculated   0.816 HSP70          4     4 large    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.318 RBOHD          4     4 small    
## 2 measurement noninoculated inoculated  -0.653 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated  -1.26  CPI8           4     4 large    
## 4 measurement noninoculated inoculated  -1.60  CAB            4     4 large    
## 5 measurement noninoculated inoculated  -0.499 BGLU2          4     4 small    
## 6 measurement noninoculated inoculated  -2.82  HSP70          4     4 large    
## 7 measurement noninoculated inoculated  -2.43  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -3.09  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result$stat.test))]
print(res[res$p.adj.signif != 'ns', ])
##   transcript        group1     group2          p      p.adj p.adj.signif
## 6      HSP70 noninoculated inoculated 0.01428571 0.03809524            *
## 7       PTI5 noninoculated inoculated 0.01428571 0.03809524            *
## 8     13-LOX noninoculated inoculated 0.01428571 0.03809524            *
cat("", file = fr, append = TRUE, sep = "\n")
output_text = "permutational t-test"
cat(output_text, file = fr, append = TRUE, sep = "\n")
header = paste(colnames(res), collapse = "\t")
cat(header, file = fr, append = TRUE, sep = "\n")
output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")


ggsave(
  filename = file.path(my_dir, "Test_roots.216_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_roots.216_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

4.3.3 shoots 218

max_per_transcript = shoots.218 %>%
  dplyr::group_by(transcript) %>%
  dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
group2_names = setdiff(max_per_transcript$transcript, group1_names)
y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
y_scales2 = build_y_scales_for(group2_names, max_per_transcript)

result = test_and_plot(data_long_raw = shoots.218
                       , myorder = myorder
                       , pal = pal
                       , what = "shoots 218"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 RBOHD             0.658               0.495             0.634  
## 2     8 PR1B              0.000437            0.000450          0.00552
## 3     8 CPI8              0.0280              0.0266            0.0614 
## 4     8 CAB               0.658               0.489             0.604  
## 5     8 BGLU2             0.522               0.489             0.604  
## 6     8 HSP70             0.658               0.505             0.604  
## 7     8 PTI5              0.522               0.489             0.604  
## 8     8 13-LOX            0.522               0.489             0.604  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6  0.255    0.631
## 2 PR1B           1     6  0.585    0.473
## 3 CPI8           1     6  3.74     0.101
## 4 CAB            1     6  0.0364   0.855
## 5 BGLU2          1     6  0.00507  0.946
## 6 HSP70          1     6  0.000737 0.979
## 7 PTI5           1     6  0.0321   0.864
## 8 13-LOX         1     6  1.91     0.216
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 RBOHD          1     6  0.255    0.631
## 2 PR1B           1     6  0.585    0.473
## 3 CPI8           1     6  3.74     0.101
## 4 CAB            1     6  0.0364   0.855
## 5 BGLU2          1     6  0.00507  0.946
## 6 HSP70          1     6  0.000737 0.979
## 7 PTI5           1     6  0.0321   0.864
## 8 13-LOX         1     6  1.91     0.216
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 RBOHD         0.715   0.398          1 Fligner-Killeen test of homogeneity of…
## 2 PR1B          0.303   0.582          1 Fligner-Killeen test of homogeneity of…
## 3 CPI8          5.03    0.0249         1 Fligner-Killeen test of homogeneity of…
## 4 CAB           0.0947  0.758          1 Fligner-Killeen test of homogeneity of…
## 5 BGLU2         0.623   0.430          1 Fligner-Killeen test of homogeneity of…
## 6 HSP70         0.0955  0.757          1 Fligner-Killeen test of homogeneity of…
## 7 PTI5          0.0955  0.757          1 Fligner-Killeen test of homogeneity of…
## 8 13-LOX        0.709   0.400          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.408 RBOHD          4     4 moderate 
## 2 measurement noninoculated inoculated   0.204 PR1B           4     4 small    
## 3 measurement noninoculated inoculated   0     CPI8           4     4 small    
## 4 measurement noninoculated inoculated   0.306 CAB            4     4 moderate 
## 5 measurement noninoculated inoculated   0.612 BGLU2          4     4 large    
## 6 measurement noninoculated inoculated   0.612 HSP70          4     4 large    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.612 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -1.23  RBOHD          4     4 large    
## 2 measurement noninoculated inoculated  -0.702 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated   0.436 CPI8           4     4 small    
## 4 measurement noninoculated inoculated  -0.615 CAB            4     4 moderate 
## 5 measurement noninoculated inoculated  -1.79  BGLU2          4     4 large    
## 6 measurement noninoculated inoculated  -1.77  HSP70          4     4 large    
## 7 measurement noninoculated inoculated  -4.62  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -1.37  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result$stat.test))]
print(res[res$p.adj.signif != 'ns', ])
## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)
cat("", file = fr, append = TRUE, sep = "\n")
output_text = "permutational t-test"
cat(output_text, file = fr, append = TRUE, sep = "\n")
header = paste(colnames(res), collapse = "\t")
cat(header, file = fr, append = TRUE, sep = "\n")
output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")



ggsave(
  filename = file.path(my_dir, "Test_shoots.218_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_shoots.218_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

4.3.4 roots 218

max_per_transcript = roots.218 %>%
  dplyr::group_by(transcript) %>%
  dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
group2_names = setdiff(max_per_transcript$transcript, group1_names)
y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
y_scales2 = build_y_scales_for(group2_names, max_per_transcript)

result = test_and_plot(data_long_raw = roots.218
                       , myorder = myorder
                       , pal = pal
                       , what = "roots 218"
                       , plot_gene_expression_func = plot_gene_expression
                       , groupvars = c("Treatment", "transcript")
                       , y_scales1
                       , y_scales2
                       )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 8 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 RBOHD              0.271               0.322            0.573  
## 2     8 PR1B               0.00339             0.00629          0.00158
## 3     8 CPI8               0.0799              0.115            0.296  
## 4     8 CAB                0.0650              0.0596           0.0195 
## 5     8 BGLU2              0.320               0.348            0.573  
## 6     8 HSP70              0.320               0.335            0.539  
## 7     8 PTI5               0.107               0.121            0.296  
## 8     8 13-LOX             0.327               0.378            0.573  
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 RBOHD          1     6    8.01   0.0300 
## 2 PR1B           1     6    0.668  0.445  
## 3 CPI8           1     6    2.55   0.161  
## 4 CAB            1     6    6.52   0.0433 
## 5 BGLU2          1     6    1.96   0.211  
## 6 HSP70          1     6    0.0221 0.887  
## 7 PTI5           1     6   24.4    0.00260
## 8 13-LOX         1     6    4.36   0.0818 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 8 × 5
##   transcript   df1   df2 statistic       p
##   <fct>      <int> <int>     <dbl>   <dbl>
## 1 RBOHD          1     6    8.01   0.0300 
## 2 PR1B           1     6    0.668  0.445  
## 3 CPI8           1     6    2.55   0.161  
## 4 CAB            1     6    6.52   0.0433 
## 5 BGLU2          1     6    1.96   0.211  
## 6 HSP70          1     6    0.0221 0.887  
## 7 PTI5           1     6   24.4    0.00260
## 8 13-LOX         1     6    4.36   0.0818 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 8 × 5
## # Groups:   transcript [8]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 RBOHD       5.03      0.0249         1 Fligner-Killeen test of homogeneity of…
## 2 PR1B        0.123     0.726          1 Fligner-Killeen test of homogeneity of…
## 3 CPI8        2.85      0.0913         1 Fligner-Killeen test of homogeneity of…
## 4 CAB         4.96      0.0259         1 Fligner-Killeen test of homogeneity of…
## 5 BGLU2       1.75      0.186          1 Fligner-Killeen test of homogeneity of…
## 6 HSP70       0.000272  0.987          1 Fligner-Killeen test of homogeneity of…
## 7 PTI5        5.00      0.0253         1 Fligner-Killeen test of homogeneity of…
## 8 13-LOX      1.76      0.185          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.306 RBOHD          4     4 moderate 
## 2 measurement noninoculated inoculated   0.408 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated   0.510 CPI8           4     4 large    
## 4 measurement noninoculated inoculated   0.816 CAB            4     4 large    
## 5 measurement noninoculated inoculated   0.408 BGLU2          4     4 moderate 
## 6 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 7 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 8 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 8 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.928 RBOHD          4     4 large    
## 2 measurement noninoculated inoculated  -0.769 PR1B           4     4 moderate 
## 3 measurement noninoculated inoculated  -1.44  CPI8           4     4 large    
## 4 measurement noninoculated inoculated  -1.94  CAB            4     4 large    
## 5 measurement noninoculated inoculated   0.906 BGLU2          4     4 large    
## 6 measurement noninoculated inoculated  -0.560 HSP70          4     4 moderate 
## 7 measurement noninoculated inoculated  -2.12  PTI5           4     4 large    
## 8 measurement noninoculated inoculated  -2.81  13-LOX         4     4 large

res = result$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result$stat.test))]
print(res[res$p.adj.signif != 'ns', ])
##   transcript        group1     group2          p      p.adj p.adj.signif
## 4        CAB noninoculated inoculated 0.01428571 0.03809524            *
## 7       PTI5 noninoculated inoculated 0.01428571 0.03809524            *
## 8     13-LOX noninoculated inoculated 0.01428571 0.03809524            *
cat("", file = fr, append = TRUE, sep = "\n")
output_text = "permutational t-test"
cat(output_text, file = fr, append = TRUE, sep = "\n")
header = paste(colnames(res), collapse = "\t")
cat(header, file = fr, append = TRUE, sep = "\n")
output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
cat(output_text, file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")

ggsave(
  filename = file.path(my_dir, "Test_roots.218_1.pdf"),
  plot = result$plot1,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 3,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)

ggsave(
  filename = file.path(my_dir, "Test_roots.218_2.pdf"),
  plot = result$plot2,
  device = pdf,
  path = NULL,
  scale = 1,
  width = 9,
  height = 8,
  units = c("in"),
  dpi = 900,
  limitsize = TRUE,
  bg = NULL
)
keep = c("magrittr", "ggplot2", "%!in%", "pal", "my_dir", "fr", "group1_names")
all_objs = ls()
# Objects that are functions or named in keep
keep = c(lsf.str(), intersect(all_objs, keep))
rm(list = setdiff(all_objs, keep))

5 Exp1

5.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "AMFcolonisation.xlsx" "qPCRmycorrhiza.xlsx"  "README.md"           
## [4] "Table S1.xlsx"        "Table S2.xlsx"        "Table S7.xlsx"
fn = 'Table S7.xlsx'

PS218 = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = 'Exp1',
                            startRow = 1,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

data.table::setDT(PS218)
PS218[, Sample.ID := NULL]
PS218[, Genotype := NULL]
PS218[, Strain := NULL]
PS218$Tissue = as.factor(trimws(PS218$Tissue))
PS218$Treatment = factor(trimws(PS218$Treatment), levels = c("noninoculated", "inoculated"))

table(PS218$Time)
## 
## 2 h 4 h 6 h 
##  16  16  16

5.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("BGLU2", "HSP70", "PTI5", "13-LOX")

5.3 results

cat("EXP1", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")


run_analysis_for_time <- function(data
                                  , time_point
                                  , myorder
                                  , pal
                                  , groupvars
                                  , measurevar
                                  , my_dir
                                  , plot_gene_expression
                                  , test_and_plot) {

  temp = data[data$Time == time_point, ]
  temp[, Time := NULL]

  shoots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$shoots
  roots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$roots
  
  max_per_transcript = shoots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  group2_names = setdiff(max_per_transcript$transcript, group1_names)
  y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales2 = build_y_scales_for(group2_names, max_per_transcript)
  max_per_transcript = roots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  y_scales3 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales4 = build_y_scales_for(group2_names, max_per_transcript)


  

  result_shoots = test_and_plot(data_long_raw = shoots,
                                myorder = myorder,
                                pal = pal,
                                what = paste("shoots", time_point),
                                plot_gene_expression_func = plot_gene_expression,
                                groupvars = groupvars,
                                y_scales1 = y_scales1,
                                y_scales2 = y_scales2)
  res = result_shoots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_shoots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Exp1_shoots.", gsub(" ", "", time_point), "_1.pdf")),
         plot = result_shoots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Exp1_shoots.", gsub(" ", "", time_point), "_2.pdf")),
         plot = result_shoots$plot2, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  

  result_roots = test_and_plot(data_long_raw = roots,
                               myorder = myorder,
                               pal = pal,
                               what = paste("roots", time_point),
                               plot_gene_expression_func = plot_gene_expression,
                               groupvars = groupvars,
                               y_scales1 = y_scales3,
                               y_scales2 = y_scales4)
  res = result_roots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_roots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  
  
  ggsave(filename = file.path(my_dir, paste0("Exp1_roots.", gsub(" ", ".", time_point), "_1.pdf")),
         plot = result_roots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Exp1_roots.", gsub(" ", ".", time_point), "_2.pdf")),
         plot = result_roots$plot2, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  

  list(shoots = result_shoots, roots = result_roots)
}


results_2h = run_analysis_for_time(data = PS218
                                   , time_point = "2 h"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.773               0.687            0.574
## 2     8 HSP70                0.773               0.687            0.574
## 3     8 PTI5                 0.773               0.687            0.631
## 4     8 13-LOX               0.773               0.696            0.689
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6    1.91   0.216
## 2 HSP70          1     6    0.0663 0.805
## 3 PTI5           1     6    0.156  0.707
## 4 13-LOX         1     6    0.0281 0.872
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6    1.91   0.216
## 2 HSP70          1     6    0.0663 0.805
## 3 PTI5           1     6    0.156  0.707
## 4 13-LOX         1     6    0.0281 0.872
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2          1.76    0.185         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70          0.619   0.431         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5           0.303   0.582         1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX         0.120   0.729         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.510 BGLU2          4     4 large    
## 2 measurement noninoculated inoculated   0.102 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.612 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0     13-LOX         4     4 small    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  1.13   BGLU2          4     4 large     
## 2 measurement noninoculated inoculated -0.176  HSP70          4     4 negligible
## 3 measurement noninoculated inoculated -1.19   PTI5           4     4 large     
## 4 measurement noninoculated inoculated  0.0435 13-LOX         4     4 negligible

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     7 BGLU2               0.252              NA                 0.717
## 2     8 HSP70               0.0417              0.0347            0.148
## 3     8 PTI5                0.0538              0.0606            0.356
## 4     8 13-LOX              0.0417              0.0347            0.144
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     2.63  0.166   
## 2 HSP70          1     6    38.9   0.000785
## 3 PTI5           1     6     4.06  0.0906  
## 4 13-LOX         1     6     0.370 0.566   
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     2.63  0.166   
## 2 HSP70          1     6    38.9   0.000785
## 3 PTI5           1     6     4.06  0.0906  
## 4 13-LOX         1     6     0.370 0.566   
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         0.923   0.337          1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         4.99    0.0255         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          5.00    0.0253         1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        0.0947  0.758          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.267 BGLU2          4     3 small    
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.789 BGLU2          4     3 moderate 
## 2 measurement noninoculated inoculated  -0.823 HSP70          4     4 large    
## 3 measurement noninoculated inoculated  -1.99  PTI5           4     4 large    
## 4 measurement noninoculated inoculated -11.7   13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 3       PTI5 noninoculated inoculated 0.01428571 0.02857143            *
## 4     13-LOX noninoculated inoculated 0.01428571 0.02857143            *

results_4h = run_analysis_for_time(data = PS218
                                   , time_point = "4 h"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2               0.424                0.492           0.703 
## 2     8 HSP70               0.424                0.492           0.723 
## 3     7 PTI5                0.0396              NA               0.0196
## 4     8 13-LOX              0.424                0.492           0.723 
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    7.26   0.0358
## 2 HSP70          1     6    2.06   0.201 
## 3 PTI5           1     5    3.62   0.115 
## 4 13-LOX         1     6    0.0744 0.794 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    7.26   0.0358
## 2 HSP70          1     6    2.06   0.201 
## 3 PTI5           1     5    3.62   0.115 
## 4 13-LOX         1     6    0.0744 0.794 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         2.84    0.0920         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         2.87    0.0903         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          0.935   0.334          1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        0.0945  0.759          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.408 BGLU2          4     4 moderate 
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.668 PTI5           4     3 large    
## 4 measurement noninoculated inoculated   0.510 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  1.14   BGLU2          4     4 large     
## 2 measurement noninoculated inoculated -0.0196 HSP70          4     4 negligible
## 3 measurement noninoculated inoculated -1.41   PTI5           4     3 large     
## 4 measurement noninoculated inoculated -0.802  13-LOX         4     4 large

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.343               0.350            0.455
## 2     8 HSP70                0.343               0.350            0.543
## 3     7 PTI5                 0.105              NA                0.252
## 4     8 13-LOX               0.343               0.350            0.453
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    0.0398 0.849 
## 2 HSP70          1     6    0.494  0.508 
## 3 PTI5           1     5   13.7    0.0140
## 4 13-LOX         1     6    2.52   0.164 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    0.0398 0.849 
## 2 HSP70          1     6    0.494  0.508 
## 3 PTI5           1     5   13.7    0.0140
## 4 13-LOX         1     6    2.52   0.164 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         0.0947  0.758          1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         0.715   0.398          1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          4.01    0.0453         1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        2.85    0.0913         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.204 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0.714 HSP70          4     4 large    
## 3 measurement noninoculated inoculated   0.802 PTI5           3     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated   0.101 BGLU2          4     4 negligible
## 2 measurement noninoculated inoculated  -2.66  HSP70          4     4 large     
## 3 measurement noninoculated inoculated  -1.62  PTI5           3     4 large     
## 4 measurement noninoculated inoculated  -2.65  13-LOX         4     4 large

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

results_6h = run_analysis_for_time(data = PS218
                                   , time_point = "6 h"
                                   , myorder
                                   , pal
                                   , groupvars
                                   , measurevar
                                   , my_dir
                                   , plot_gene_expression
                                   , test_and_plot)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.410               0.420            0.585
## 2     8 HSP70                0.932               0.825            0.780
## 3     8 PTI5                 0.479               0.696            0.780
## 4     8 13-LOX               0.932               0.825            0.780
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6     9.33  0.0224
## 2 HSP70          1     6     1.76  0.233 
## 3 PTI5           1     6     0.632 0.457 
## 4 13-LOX         1     6     0.157 0.705 
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6     9.33  0.0224
## 2 HSP70          1     6     1.76  0.233 
## 3 PTI5           1     6     0.632 0.457 
## 4 13-LOX         1     6     0.157 0.705 
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         4.99    0.0255         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         1.47    0.226          1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          1.76    0.185          1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        0.0937  0.760          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.204 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.714 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.510 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  -0.783 BGLU2          4     4 moderate  
## 2 measurement noninoculated inoculated  -0.167 HSP70          4     4 negligible
## 3 measurement noninoculated inoculated  -2.26  PTI5           4     4 large     
## 4 measurement noninoculated inoculated  -1.01  13-LOX         4     4 large

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.163               0.184            0.127
## 2     8 HSP70                0.514               0.450            0.330
## 3     8 PTI5                 0.163               0.184            0.127
## 4     8 13-LOX               0.163               0.185            0.330
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    1.48   0.269 
## 2 HSP70          1     6    0.0375 0.853 
## 3 PTI5           1     6    1.80   0.228 
## 4 13-LOX         1     6    6.91   0.0391
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic      p
##   <fct>      <int> <int>     <dbl>  <dbl>
## 1 BGLU2          1     6    1.48   0.269 
## 2 HSP70          1     6    0.0375 0.853 
## 3 PTI5           1     6    1.80   0.228 
## 4 13-LOX         1     6    6.91   0.0391
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         0.123   0.726          1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         0.0947  0.758          1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          1.47    0.226          1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        5.00    0.0253         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.714 BGLU2          4     4 large    
## 2 measurement noninoculated inoculated   0.306 HSP70          4     4 moderate 
## 3 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -2.10  BGLU2          4     4 large    
## 2 measurement noninoculated inoculated  -0.886 HSP70          4     4 large    
## 3 measurement noninoculated inoculated  -4.39  PTI5           4     4 large    
## 4 measurement noninoculated inoculated  -2.64  13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 3       PTI5 noninoculated inoculated 0.01428571 0.02857143            *
## 4     13-LOX noninoculated inoculated 0.01428571 0.02857143            *

keep = c("magrittr", "ggplot2", "%!in%", "pal", "my_dir", "fr", "group1_names")
all_objs = ls()
# Objects that are functions or named in keep
keep = c(lsf.str(), intersect(all_objs, keep))
rm(list = setdiff(all_objs, keep))

6 Exp2

6.1 input

fp = file.path('..', "input")
list.files(fp)
## [1] "AMFcolonisation.xlsx" "qPCRmycorrhiza.xlsx"  "README.md"           
## [4] "Table S1.xlsx"        "Table S2.xlsx"        "Table S7.xlsx"
fn = 'Table S7.xlsx'

PS218 = openxlsx::read.xlsx(xlsxFile = file.path(fp, fn),
                            sheet = 'Exp2',
                            startRow = 1,
                            colNames = TRUE,
                            rowNames = FALSE,
                            detectDates = FALSE,
                            skipEmptyRows = TRUE,
                            skipEmptyCols = TRUE,
                            rows = NULL,
                            cols = NULL,
                            check.names = FALSE,
                            sep.names = ".",
                            namedRegion = NULL,
                            na.strings = "NA",
                            fillMergedCells = FALSE)

data.table::setDT(PS218)
PS218[, Sample.ID := NULL]
PS218[, Genotype := NULL]
PS218[, Strain := NULL]
PS218$Tissue = as.factor(trimws(PS218$Tissue))
PS218$Treatment = factor(trimws(PS218$Treatment), levels = c("noninoculated", "inoculated"))

table(PS218$Time)
## 
##    1 h  1 min 15 min    2 h 30 min 
##      8      8      8      8      8

6.2 params

groupvars = c("Tissue", "Treatment", "transcript")
measurevar = 'measurement'
dodge = position_dodge(width = 0.5)

myorder = c("BGLU2", "HSP70", "PTI5", "13-LOX")

cat("EXP2", file = fr, append = TRUE, sep = "\n")
cat("", file = fr, append = TRUE, sep = "\n")

6.3 results

run_analysis_for_time <- function(data
                                  , time_point
                                  , myorder
                                  , pal
                                  , groupvars
                                  , measurevar
                                  , y_scales1
                                  , y_scales2
                                  , my_dir
                                  , plot_gene_expression
                                  , test_and_plot
                                  ) {

  temp = data[data$Time == time_point, ]
  temp[, Time := NULL]

  roots = process_data(temp, groupvars, measurevar, scale_treatment = "noninoculated")$roots
  
  max_per_transcript = roots %>%
    dplyr::group_by(transcript) %>%
    dplyr::summarise(max_scaled = max(scaled, na.rm = TRUE), .groups = "drop")
  group2_names = setdiff(max_per_transcript$transcript, group1_names)
  y_scales1 = build_y_scales_for(group1_names, max_per_transcript)
  y_scales2 = build_y_scales_for(group2_names, max_per_transcript)
  
  

  result_roots = test_and_plot(data_long_raw = roots,
                               myorder = myorder,
                               pal = pal,
                               what = paste("roots", time_point),
                               plot_gene_expression_func = plot_gene_expression,
                               groupvars = groupvars,
                               y_scales1,
                               y_scales2)
  res = result_roots$stat.test[, grep("transcript|group1|group2|^p$|p\\.", colnames(result_roots$stat.test))]
  print(res[res$p.adj.signif != 'ns', ])
  
  cat("", file = fr, append = TRUE, sep = "\n")
  output_text = "permutational t-test"
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  header = paste(colnames(res), collapse = "\t")
  cat(header, file = fr, append = TRUE, sep = "\n")
  output_text = apply(as.data.frame(tibble::as_tibble(res)), 1, function(row) paste(row, collapse = "\t"))
  cat(output_text, file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")
  cat("", file = fr, append = TRUE, sep = "\n")


  
  ggsave(filename = file.path(my_dir, paste0("Exp2_roots.", gsub(" ", ".", time_point), "_1.pdf")),
         plot = result_roots$plot1, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  ggsave(filename = file.path(my_dir, paste0("Exp2_roots.", gsub(" ", ".", time_point), "_2.pdf")),
         plot = result_roots$plot2, device = pdf, width = 3, height = 8, units = "in", dpi = 900)
  

  list(roots = result_roots)
}



results_1min = run_analysis_for_time(data = PS218
                                     , time_point = "1 min"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.556               0.328            0.343
## 2     8 HSP70                0.774               0.511            0.400
## 3     8 PTI5                 0.575               0.328            0.400
## 4     8 13-LOX               0.556               0.328            0.400
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6  1.63     0.249
## 2 HSP70          1     6  0.530    0.494
## 3 PTI5           1     6  0.000531 0.982
## 4 13-LOX         1     6  0.652    0.450
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6  1.63     0.249
## 2 HSP70          1     6  0.530    0.494
## 3 PTI5           1     6  0.000531 0.982
## 4 13-LOX         1     6  0.652    0.450
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         0.123    0.726         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         0.710    0.399         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          0.0945   0.759         1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        0.304    0.581         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.612 BGLU2          4     4 large    
## 2 measurement noninoculated inoculated   0.204 HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.102 PTI5           4     4 small    
## 4 measurement noninoculated inoculated   0.204 13-LOX         4     4 small    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated  -1.05  BGLU2          4     4 large     
## 2 measurement noninoculated inoculated  -0.798 HSP70          4     4 moderate  
## 3 measurement noninoculated inoculated   0.430 PTI5           4     4 small     
## 4 measurement noninoculated inoculated   0.137 13-LOX         4     4 negligible

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

results_15min = run_analysis_for_time(data = PS218
                                     , time_point = "15 min"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.480               0.666            0.948
## 2     8 HSP70                0.920               0.910            0.948
## 3     8 PTI5                 0.920               0.910            0.948
## 4     8 13-LOX               0.920               0.910            0.948
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     2.83  0.143
## 2 HSP70          1     6     0.795 0.407
## 3 PTI5           1     6     0.890 0.382
## 4 13-LOX         1     6     3.67  0.104
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     2.83  0.143
## 2 HSP70          1     6     0.795 0.407
## 3 PTI5           1     6     0.890 0.382
## 4 13-LOX         1     6     3.67  0.104
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2           2.84  0.0920         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70           1.47  0.226          1 Fligner-Killeen test of homogeneity of…
## 3 PTI5            1.46  0.226          1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX          1.46  0.227          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.306 BGLU2          4     4 moderate 
## 2 measurement noninoculated inoculated   0.816 HSP70          4     4 large    
## 3 measurement noninoculated inoculated   0.408 PTI5           4     4 moderate 
## 4 measurement noninoculated inoculated   0.204 13-LOX         4     4 small    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.662 BGLU2          4     4 moderate 
## 2 measurement noninoculated inoculated   2.07  HSP70          4     4 large    
## 3 measurement noninoculated inoculated  -0.883 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.623 13-LOX         4     4 moderate

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

results_30min = run_analysis_for_time(data = PS218
                                     , time_point = "30 min"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2                0.734               0.672            0.579
## 2     8 HSP70                0.734               0.672            0.579
## 3     8 PTI5                 0.734               0.672            0.579
## 4     8 13-LOX               0.937               0.921            0.788
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.752 0.419
## 2 HSP70          1     6     1.64  0.248
## 3 PTI5           1     6     0.310 0.598
## 4 13-LOX         1     6     1.36  0.288
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.752 0.419
## 2 HSP70          1     6     1.64  0.248
## 3 PTI5           1     6     0.310 0.598
## 4 13-LOX         1     6     1.36  0.288
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2          1.75    0.186         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70          0.715   0.398         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5           0.306   0.580         1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX         0.619   0.432         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.102 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0.408 HSP70          4     4 moderate 
## 3 measurement noninoculated inoculated   0.612 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.204 13-LOX         4     4 small    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude 
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>     
## 1 measurement noninoculated inoculated -0.0596 BGLU2          4     4 negligible
## 2 measurement noninoculated inoculated -0.800  HSP70          4     4 moderate  
## 3 measurement noninoculated inoculated -1.16   PTI5           4     4 large     
## 4 measurement noninoculated inoculated -0.497  13-LOX         4     4 small

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

results_1h = run_analysis_for_time(data = PS218
                                     , time_point = "1 h"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     7 BGLU2               0.363                   NA            0.718
## 2     7 HSP70               0.0973                  NA            0.188
## 3     7 PTI5                0.0782                  NA            0.188
## 4     7 13-LOX              0.363                   NA            0.718
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     93.6  0.000200
## 2 HSP70          1     5      5.10 0.0736  
## 3 PTI5           1     5      1.07 0.348   
## 4 13-LOX         1     5      1.85 0.232   
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic        p
##   <fct>      <int> <int>     <dbl>    <dbl>
## 1 BGLU2          1     5     93.6  0.000200
## 2 HSP70          1     5      5.10 0.0736  
## 3 PTI5           1     5      1.07 0.348   
## 4 13-LOX         1     5      1.85 0.232   
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2          4.01   0.0453         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70          4.05   0.0441         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5           0.662  0.416          1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX         1.89   0.170          1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0     BGLU2          3     4 small    
## 2 measurement noninoculated inoculated   0.401 HSP70          3     4 moderate 
## 3 measurement noninoculated inoculated   0.668 PTI5           3     4 large    
## 4 measurement noninoculated inoculated   0.802 13-LOX         3     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated  -0.393 BGLU2          3     4 small    
## 2 measurement noninoculated inoculated  -0.912 HSP70          3     4 large    
## 3 measurement noninoculated inoculated  -1.35  PTI5           3     4 large    
## 4 measurement noninoculated inoculated  -3.09  13-LOX         3     4 large

## [1] transcript   group1       group2       p            p.adj       
## [6] p.adj.signif
## <0 rows> (or 0-length row.names)

results_2h = run_analysis_for_time(data = PS218
                                     , time_point = "2 h"
                                     , myorder
                                     , pal
                                     , groupvars
                                     , measurevar
                                     , y_scales1
                                     , y_scales2
                                     , my_dir
                                     , plot_gene_expression
                                     , test_and_plot
                                     )

## ####  ####  
## Distribution tests
## ####  ####  
## # A tibble: 4 × 7
##       n transcript Shapiro_Wilk_BH Anderson_Darling_BH Lilliefors_KS_BH
##   <int> <fct>                <dbl>               <dbl>            <dbl>
## 1     8 BGLU2               0.215               0.290            0.755 
## 2     8 HSP70               0.0287              0.0275           0.0291
## 3     8 PTI5                0.0959              0.127            0.400 
## 4     8 13-LOX              0.0287              0.0275           0.0914
## # ℹ 2 more variables: Jarque_Bera_BH <dbl>, DAgostino_Skewness_BH <dbl>
## ####  ####  
## Quantile-Quantile plots
## ####  ####

## ####  ####  
## Test for homogeneity of variance across groups
## ####  ####  
## ####  ####  
## Levene
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.151 0.711
## 2 HSP70          1     6     0.926 0.373
## 3 PTI5           1     6     2.14  0.194
## 4 13-LOX         1     6     0.634 0.456
## ####  ####  
## Brown-Forsythe
## ####  ####  
## # A tibble: 4 × 5
##   transcript   df1   df2 statistic     p
##   <fct>      <int> <int>     <dbl> <dbl>
## 1 BGLU2          1     6     0.151 0.711
## 2 HSP70          1     6     0.926 0.373
## 3 PTI5           1     6     2.14  0.194
## 4 13-LOX         1     6     0.634 0.456
## ####  ####  
## Fligner
## ####  ####  
## # A tibble: 4 × 5
## # Groups:   transcript [4]
##   transcript statistic p.value parameter method                                 
##   <fct>          <dbl>   <dbl>     <dbl> <chr>                                  
## 1 BGLU2         0.0945   0.759         1 Fligner-Killeen test of homogeneity of…
## 2 HSP70         0.120    0.729         1 Fligner-Killeen test of homogeneity of…
## 3 PTI5          1.45     0.228         1 Fligner-Killeen test of homogeneity of…
## 4 13-LOX        1.46     0.226         1 Fligner-Killeen test of homogeneity of…
## ####  ####  
## Wilcoxon effect size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.102 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated   0     HSP70          4     4 small    
## 3 measurement noninoculated inoculated   0.816 PTI5           4     4 large    
## 4 measurement noninoculated inoculated   0.816 13-LOX         4     4 large    
## ####  ####  
## Cohen's d Measure of Effect Size
## ####  ####  
## # A tibble: 4 × 8
##   .y.         group1        group2     effsize transcript    n1    n2 magnitude
## * <chr>       <chr>         <chr>        <dbl> <fct>      <int> <int> <ord>    
## 1 measurement noninoculated inoculated   0.413 BGLU2          4     4 small    
## 2 measurement noninoculated inoculated  -0.442 HSP70          4     4 small    
## 3 measurement noninoculated inoculated  -2.37  PTI5           4     4 large    
## 4 measurement noninoculated inoculated -12.6   13-LOX         4     4 large

##   transcript        group1     group2          p      p.adj p.adj.signif
## 3       PTI5 noninoculated inoculated 0.01428571 0.02857143            *
## 4     13-LOX noninoculated inoculated 0.01428571 0.02857143            *

sessionInfo()
## R version 4.4.1 (2024-06-14 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8 
## [2] LC_CTYPE=English_United Kingdom.utf8   
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.utf8    
## 
## time zone: Europe/Ljubljana
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_4.0.0  magrittr_2.0.3
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.3             Rdpack_2.6.4          sandwich_3.1-1       
##   [4] MKinfer_1.2           rlang_1.1.5           multcomp_1.4-28      
##   [7] miceadds_3.18-36      tseries_0.10-58       matrixStats_1.5.0    
##  [10] compiler_4.4.1        vctrs_0.6.5           quadprog_1.5-8       
##  [13] pkgconfig_2.0.3       shape_1.4.6.1         crayon_1.5.3         
##  [16] fastmap_1.2.0         backports_1.5.0       labeling_0.4.3       
##  [19] utf8_1.2.6            rmarkdown_2.30        exactRankTests_0.8-35
##  [22] nloptr_2.2.1          purrr_1.0.4           xfun_0.53            
##  [25] glmnet_4.1-10         modeltools_0.2-24     jomo_2.7-6           
##  [28] cachem_1.1.0          jsonlite_2.0.0        gmp_0.7-5            
##  [31] pan_1.9               broom_1.0.10          parallel_4.4.1       
##  [34] R6_2.6.1              coin_1.4-3            bslib_0.9.0          
##  [37] stringi_1.8.7         RColorBrewer_1.1-3    rpart_4.1.24         
##  [40] car_3.1-3             boot_1.3-31           jquerylib_0.1.4      
##  [43] Rcpp_1.0.14           iterators_1.0.14      knitr_1.50           
##  [46] zoo_1.8-14            nnet_7.3-20           Matrix_1.7-1         
##  [49] splines_4.4.1         tidyselect_1.2.1      rstudioapi_0.17.1    
##  [52] dichromat_2.0-0.1     abind_1.4-8           yaml_2.3.10          
##  [55] codetools_0.2-20      curl_6.2.2            arrangements_1.1.9   
##  [58] lattice_0.22-6        tibble_3.3.0          quantmod_0.4.28      
##  [61] withr_3.0.2           S7_0.2.0              evaluate_1.0.5       
##  [64] moments_0.14.1        survival_3.8-3        zip_2.3.3            
##  [67] xts_0.14.1            pillar_1.11.1         ggpubr_0.6.1         
##  [70] carData_3.0-5         mice_3.18.0           nortest_1.0-4        
##  [73] foreach_1.5.2         stats4_4.4.1          reformulas_0.4.1     
##  [76] generics_0.1.4        TTR_0.24.4            scales_1.4.0         
##  [79] minqa_1.2.8           glue_1.8.0            tools_4.4.1          
##  [82] data.table_1.17.8     MKdescr_0.9           lme4_1.1-37          
##  [85] openxlsx_4.2.8        ggsignif_0.6.4        mvtnorm_1.3-3        
##  [88] cowplot_1.2.0         grid_4.4.1            mitools_2.4          
##  [91] tidyr_1.3.1           rbibutils_2.3         libcoin_1.0-10       
##  [94] nlme_3.1-166          Formula_1.2-5         cli_3.6.3            
##  [97] dplyr_1.1.4           ggh4x_0.3.1           gtable_0.3.6         
## [100] rstatix_0.7.2         sass_0.4.10           digest_0.6.37        
## [103] TH.data_1.1-4         farver_2.1.2          htmltools_0.5.8.1    
## [106] lifecycle_1.0.4       mitml_0.4-5           MASS_7.3-64